stanEx4
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # not see parameters str1..., str2,... using regex as exc='[str1,str2,...]'
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}explanatory variable obs.x have noise
xN(x0.sx),yN(b0+b1*x0,s)
ex8-0.stan
normal regression without explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}ex9.stan
normal regression with explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x10;
}
parameters{
real b0;
real b1;
real<lower=0> s;
real<lower=0> sx;
vector[N] x0;
}
model{
for(i in 1:N){
x[i]~normal(x0[i],sx);
y[i]~normal(b0+b1*x0[i],s);
}
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x0[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] x1;
vector[N1] y1;
for(i in 1:N1){
x1[i]=normal_rng(x10[i],sx);
m1[i]=b0+b1*x10[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)
n1=10
#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -33279.9
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 25 -29.4414 0.00414791 0.000567954 1 1 53
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -29.44
## b0 6.69
## b1 1.80
## s 2.64
## m0[1] 47.65
## m0[2] 8.21
## m0[3] 26.38
## m0[4] 27.20
## m0[5] 29.29
## m0[6] 22.88
## m0[7] 21.40
## m0[8] 41.24
## m0[9] 18.12
## m0[10] 16.17
## m0[11] 30.09
## m0[12] 20.15
## m0[13] 28.04
## m0[14] 36.42
## m0[15] 42.89
## m0[16] 42.80
## m0[17] 11.65
## m0[18] 8.08
## m0[19] 9.86
## m0[20] 10.35
## y0[1] 43.65
## y0[2] 4.97
## y0[3] 26.89
## y0[4] 26.51
## y0[5] 32.87
## y0[6] 24.22
## y0[7] 23.61
## y0[8] 47.38
## y0[9] 17.91
## y0[10] 15.97
## y0[11] 30.39
## y0[12] 18.75
## y0[13] 29.75
## y0[14] 36.71
## y0[15] 44.16
## y0[16] 46.78
## y0[17] 16.28
## y0[18] 8.91
## y0[19] 4.77
## y0[20] 6.08
## m1[1] 8.08
## m1[2] 12.47
## m1[3] 16.87
## m1[4] 21.27
## m1[5] 25.66
## m1[6] 30.06
## m1[7] 34.46
## m1[8] 38.86
## m1[9] 43.25
## m1[10] 47.65
## y1[1] 5.09
## y1[2] 8.37
## y1[3] 15.28
## y1[4] 20.41
## y1[5] 24.67
## y1[6] 28.37
## y1[7] 34.87
## y1[8] 43.52
## y1[9] 44.01
## y1[10] 47.75
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -30.13 -29.75 1.40 1.12 -32.99 -28.65 1.01 604 841
## b0 6.63 6.59 1.22 1.20 4.70 8.63 1.01 266 392
## b1 1.81 1.80 0.10 0.09 1.64 1.97 1.01 360 680
## s 3.01 2.95 0.58 0.50 2.23 4.10 1.00 891 986
## m0[1] 47.74 47.72 1.47 1.39 45.40 50.18 1.00 1059 1286
## m0[2] 8.16 8.11 1.15 1.12 6.35 10.06 1.01 268 391
## m0[3] 26.39 26.38 0.72 0.66 25.22 27.59 1.00 1061 875
## m0[4] 27.21 27.20 0.73 0.67 26.04 28.42 1.00 1294 838
## m0[5] 29.31 29.30 0.76 0.72 28.11 30.56 1.00 1835 1063
## m0[6] 22.88 22.84 0.72 0.69 21.76 24.09 1.01 556 730
## m0[7] 21.39 21.37 0.73 0.71 20.25 22.62 1.01 455 642
## m0[8] 41.31 41.30 1.17 1.08 39.41 43.24 1.00 1587 1208
## m0[9] 18.10 18.07 0.80 0.76 16.85 19.44 1.01 348 465
## m0[10] 16.15 16.13 0.85 0.80 14.80 17.61 1.01 316 503
## m0[11] 30.11 30.10 0.78 0.74 28.88 31.38 1.00 1996 954
## m0[12] 20.14 20.12 0.75 0.72 18.95 21.37 1.01 403 563
## m0[13] 28.06 28.05 0.74 0.68 26.87 29.27 1.00 1523 961
## m0[14] 36.47 36.47 0.97 0.92 34.88 38.05 1.00 2127 1240
## m0[15] 42.96 42.95 1.24 1.16 40.95 45.03 1.00 1414 1180
## m0[16] 42.87 42.86 1.24 1.15 40.87 44.93 1.00 1424 1164
## m0[17] 11.61 11.57 1.01 0.95 10.02 13.29 1.01 280 408
## m0[18] 8.02 7.98 1.16 1.13 6.20 9.93 1.01 268 391
## m0[19] 9.81 9.76 1.08 1.03 8.14 11.60 1.01 273 391
## m0[20] 10.31 10.26 1.06 1.01 8.65 12.07 1.01 274 395
## y0[1] 47.80 47.80 3.43 3.27 42.23 53.30 1.00 1614 1877
## y0[2] 8.15 8.10 3.27 3.06 2.96 13.87 1.00 1126 1624
## y0[3] 26.35 26.37 3.17 2.98 21.21 31.41 1.00 2127 2100
## y0[4] 27.19 27.22 3.15 2.90 22.07 32.56 1.00 1861 1963
## y0[5] 29.35 29.29 3.15 2.97 24.28 34.60 1.00 1769 1764
## y0[6] 22.86 22.89 3.14 3.05 17.69 27.78 1.00 1907 1891
## y0[7] 21.36 21.35 3.26 3.00 16.24 26.75 1.00 1782 1896
## y0[8] 41.31 41.31 3.35 3.04 35.87 46.69 1.00 1913 1971
## y0[9] 18.17 18.12 3.21 3.06 12.99 23.38 1.00 1789 1971
## y0[10] 16.30 16.37 3.16 3.12 11.05 21.33 1.00 1611 1698
## y0[11] 30.30 30.26 3.19 3.06 25.09 35.53 1.00 1956 1749
## y0[12] 20.11 20.10 3.13 3.00 14.90 25.10 1.00 1680 1956
## y0[13] 28.14 28.21 3.16 2.94 22.87 33.31 1.00 1728 1941
## y0[14] 36.57 36.56 3.26 3.17 31.30 41.83 1.00 1937 1972
## y0[15] 42.99 42.96 3.24 3.13 37.73 48.38 1.00 1778 1821
## y0[16] 42.90 42.95 3.27 3.15 37.41 48.19 1.00 2021 1700
## y0[17] 11.61 11.67 3.23 3.16 6.12 16.82 1.00 1653 1681
## y0[18] 8.10 7.97 3.28 3.06 2.74 13.50 1.00 1352 1564
## y0[19] 9.75 9.82 3.16 2.94 4.43 14.81 1.00 1345 1681
## y0[20] 10.30 10.19 3.24 3.06 5.14 15.78 1.01 973 1672
## m1[1] 8.02 7.98 1.16 1.13 6.20 9.93 1.01 268 391
## m1[2] 12.44 12.39 0.98 0.92 10.89 14.08 1.01 285 423
## m1[3] 16.85 16.82 0.83 0.79 15.54 18.26 1.01 325 493
## m1[4] 21.26 21.24 0.73 0.70 20.12 22.49 1.01 449 615
## m1[5] 25.68 25.66 0.72 0.66 24.53 26.89 1.00 912 829
## m1[6] 30.09 30.08 0.78 0.74 28.85 31.36 1.00 1991 954
## m1[7] 34.50 34.50 0.90 0.86 33.04 35.97 1.00 2270 1207
## m1[8] 38.91 38.92 1.07 0.99 37.16 40.67 1.00 1862 1190
## m1[9] 43.33 43.31 1.26 1.18 41.29 45.43 1.00 1379 1210
## m1[10] 47.74 47.72 1.47 1.39 45.40 50.18 1.00 1059 1286
## y1[1] 8.15 8.15 3.34 3.18 2.75 13.56 1.00 947 1287
## y1[2] 12.38 12.31 3.11 3.04 7.37 17.49 1.00 1367 1644
## y1[3] 16.78 16.77 3.11 3.03 11.71 22.06 1.00 1503 1785
## y1[4] 21.14 21.13 3.20 3.01 15.89 26.19 1.00 2127 1872
## y1[5] 25.77 25.86 3.11 3.00 20.75 30.78 1.00 1802 1847
## y1[6] 30.21 30.20 3.15 3.02 25.10 35.38 1.00 2052 1919
## y1[7] 34.61 34.65 3.22 3.14 29.25 39.92 1.00 1685 1692
## y1[8] 38.86 38.86 3.30 3.15 33.49 44.09 1.00 1967 1895
## y1[9] 43.36 43.33 3.31 3.15 37.97 49.02 1.00 2010 1946
## y1[10] 47.72 47.69 3.35 3.23 42.43 53.44 1.00 1611 1497
#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)
mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 4 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 4 finished in 0.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.seeMCMC(mcmc,'[m,x]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -27.18 -29.35 10.87 10.30 -42.65 -9.01 1.04 71 107
## b0 7.07 6.88 1.56 1.31 4.80 10.37 1.11 23 24
## b1 1.77 1.78 0.13 0.11 1.50 1.96 1.12 23 8
## s 1.91 1.96 1.02 1.21 0.33 3.51 1.11 26 10
## sx 1.17 1.19 0.60 0.65 0.20 2.12 1.09 31 93
## x0[1] 21.89 22.04 1.01 1.08 20.10 23.26 1.02 157 603
## x0[2] -0.15 0.11 1.25 1.05 -3.12 1.39 1.11 24 12
## x0[3] 10.15 10.18 0.84 1.02 8.93 11.45 1.05 67 1545
## x0[4] 11.51 11.48 0.69 0.53 10.39 12.71 1.03 1714 1673
## x0[5] 12.54 12.54 0.66 0.54 11.41 13.63 1.01 1087 1902
## x0[6] 9.17 9.11 0.69 0.61 8.13 10.36 1.02 145 1927
## x0[7] 6.31 6.40 1.53 1.86 3.50 8.36 1.11 25 9
## x0[8] 19.38 19.29 0.82 0.69 18.07 20.84 1.07 34 88
## x0[9] 6.10 6.15 0.77 0.68 4.83 7.34 1.06 47 275
## x0[10] 4.02 4.22 1.24 1.30 1.55 5.63 1.11 25 9
## x0[11] 12.32 12.32 0.78 0.89 11.12 13.53 1.02 175 2177
## x0[12] 8.20 8.22 0.80 0.84 6.96 9.50 1.03 306 1433
## x0[13] 11.60 11.61 0.68 0.60 10.52 12.71 1.01 514 2123
## x0[14] 17.58 17.47 1.08 1.16 16.08 19.61 1.10 27 34
## x0[15] 21.00 20.79 1.13 1.00 19.55 23.58 1.11 24 83
## x0[16] 21.33 21.11 1.30 1.33 19.67 23.94 1.12 22 35
## x0[17] 3.33 3.27 0.84 0.93 2.08 4.71 1.05 64 338
## x0[18] 0.93 0.94 0.89 0.78 -0.66 2.34 1.10 27 31
## x0[19] 1.91 1.94 0.86 0.73 0.25 3.29 1.08 30 13
## x0[20] 2.94 2.87 0.96 1.07 1.61 4.53 1.03 138 1395
## m0[1] 45.74 45.48 1.93 2.08 43.18 49.18 1.03 149 1014
## m0[2] 6.91 6.68 1.60 1.56 4.68 9.64 1.02 204 1656
## m0[3] 25.05 24.90 1.46 1.54 22.90 27.47 1.02 179 1805
## m0[4] 27.43 27.48 1.26 1.05 25.32 29.49 1.04 3948 1821
## m0[5] 29.23 29.23 1.22 1.02 27.19 31.27 1.02 3530 1961
## m0[6] 23.31 23.35 1.24 1.02 21.19 25.27 1.03 795 1743
## m0[7] 18.33 18.15 2.38 3.20 15.08 22.07 1.04 80 861
## m0[8] 41.29 41.28 1.35 1.14 39.08 43.51 1.02 2963 2674
## m0[9] 17.90 17.91 1.25 1.01 15.91 19.99 1.04 2510 1933
## m0[10] 14.27 14.09 1.76 2.12 11.87 17.13 1.03 131 840
## m0[11] 28.87 28.74 1.40 1.43 26.86 31.28 1.02 204 2092
## m0[12] 21.58 21.70 1.49 1.69 19.13 23.71 1.03 144 2182
## m0[13] 27.59 27.48 1.23 1.08 25.63 29.69 1.04 1142 1670
## m0[14] 38.10 38.28 1.66 1.80 35.16 40.34 1.03 135 899
## m0[15] 44.12 44.35 1.59 1.51 41.27 46.32 1.02 211 896
## m0[16] 44.69 44.98 1.85 1.91 41.36 47.18 1.02 170 789
## m0[17] 13.00 13.15 1.51 1.60 10.47 15.11 1.02 171 1659
## m0[18] 8.78 8.94 1.42 1.21 6.27 10.88 1.01 679 1943
## m0[19] 10.51 10.66 1.36 1.19 8.24 12.64 1.03 813 1773
## m0[20] 12.31 12.54 1.78 1.95 9.31 14.78 1.03 139 1147
## y0[1] 45.76 45.12 2.89 2.35 41.99 51.17 1.01 460 1591
## y0[2] 6.95 6.47 2.70 2.09 3.04 11.95 1.01 767 2474
## y0[3] 25.05 24.63 2.52 2.03 21.39 29.63 1.02 753 1704
## y0[4] 27.46 27.48 2.49 1.85 23.39 31.65 1.04 4033 2651
## y0[5] 29.21 29.19 2.51 1.90 25.13 33.44 1.04 3788 1906
## y0[6] 23.36 23.44 2.51 1.92 19.12 27.32 1.04 3393 2606
## y0[7] 18.33 17.62 3.26 3.10 14.32 24.31 1.03 135 1810
## y0[8] 41.35 41.27 2.56 1.97 37.14 45.72 1.03 3763 3186
## y0[9] 17.86 17.83 2.57 1.87 13.71 22.34 1.03 3113 2782
## y0[10] 14.30 13.69 2.82 2.34 10.59 19.48 1.01 364 1498
## y0[11] 28.87 28.42 2.62 2.00 25.03 33.81 1.02 928 2099
## y0[12] 21.59 22.03 2.62 2.14 16.87 25.26 1.02 679 2011
## y0[13] 27.56 27.43 2.46 1.85 23.48 31.88 1.04 3814 3358
## y0[14] 38.18 38.70 2.70 2.21 33.16 41.93 1.01 416 2176
## y0[15] 44.12 44.56 2.66 2.10 39.33 47.93 1.02 897 2533
## y0[16] 44.71 45.33 2.85 2.23 39.38 48.62 1.01 433 1892
## y0[17] 12.98 13.39 2.67 2.13 8.19 16.89 1.01 593 2684
## y0[18] 8.71 9.01 2.60 1.99 4.05 12.67 1.02 2333 2817
## y0[19] 10.54 10.74 2.58 1.99 5.97 14.66 1.05 3089 2202
## y0[20] 12.28 12.90 2.80 2.30 7.05 15.97 1.02 386 2314
## m1[1] 8.44 8.25 1.47 1.23 6.27 11.51 1.11 24 24
## m1[2] 12.76 12.62 1.19 1.03 10.95 15.16 1.11 24 26
## m1[3] 17.08 16.98 0.94 0.85 15.62 18.84 1.10 26 33
## m1[4] 21.40 21.36 0.75 0.75 20.20 22.57 1.07 39 299
## m1[5] 25.73 25.75 0.66 0.63 24.61 26.73 1.02 336 2311
## m1[6] 30.05 30.00 0.71 0.66 28.90 31.20 1.01 664 2036
## m1[7] 34.37 34.37 0.89 0.90 33.03 35.80 1.04 63 1123
## m1[8] 38.70 38.73 1.13 1.08 36.90 40.49 1.07 35 125
## m1[9] 43.02 43.07 1.40 1.29 40.68 45.18 1.09 29 50
## m1[10] 47.34 47.39 1.68 1.52 44.38 49.91 1.10 27 39
## x1[1] 0.74 0.75 1.33 0.98 -1.48 2.97 1.02 4338 1139
## x1[2] 3.24 3.23 1.31 0.98 1.05 5.45 1.02 4268 1673
## x1[3] 5.65 5.65 1.31 0.97 3.47 7.81 1.02 4201 1602
## x1[4] 8.10 8.10 1.31 0.96 5.99 10.26 1.02 4137 1727
## x1[5] 10.52 10.54 1.31 0.97 8.23 12.71 1.02 4032 1117
## x1[6] 13.00 13.00 1.32 0.95 10.81 15.23 1.02 3733 1580
## x1[7] 15.43 15.43 1.29 0.93 13.22 17.55 1.02 4008 1294
## x1[8] 17.87 17.88 1.28 0.97 15.77 19.93 1.02 3976 2072
## x1[9] 20.31 20.32 1.28 0.92 18.16 22.42 1.02 3946 1624
## x1[10] 22.77 22.78 1.34 0.96 20.48 24.98 1.02 3972 1511
## y1[1] 8.41 8.41 2.68 2.40 4.01 12.28 1.05 55 2207
## y1[2] 12.85 12.84 2.45 2.31 8.86 16.42 1.04 88 2845
## y1[3] 17.08 17.14 2.39 2.11 13.17 20.80 1.03 115 2504
## y1[4] 21.40 21.49 2.33 1.81 17.48 25.14 1.01 1359 2477
## y1[5] 25.71 25.82 2.28 1.69 21.78 29.47 1.03 2863 2433
## y1[6] 30.11 30.04 2.22 1.64 26.46 33.86 1.04 3095 2120
## y1[7] 34.38 34.29 2.35 1.85 30.64 38.51 1.01 1188 2880
## y1[8] 38.71 38.65 2.39 2.10 35.02 42.74 1.02 149 2805
## y1[9] 42.98 42.98 2.56 2.36 39.12 47.33 1.03 102 2725
## y1[10] 47.37 47.38 2.75 2.54 43.23 51.84 1.04 70 2183
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
outlier
ex10.stan
objective variable have outlier, y~cauchy(b0+b1*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~cauchy(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=cauchy_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=cauchy_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -29125.9
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 28 -33.6454 0.000773142 0.000329433 1 1 59
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -33.65
## b0 6.46
## b1 1.85
## s 3.26
## m0[1] 12.01
## m0[2] 14.93
## m0[3] 15.25
## m0[4] 21.16
## m0[5] 16.82
## m0[6] 18.91
## m0[7] 14.07
## m0[8] 18.98
## m0[9] 10.35
## m0[10] 13.62
## m0[11] 11.20
## m0[12] 13.32
## m0[13] 23.08
## m0[14] 8.65
## m0[15] 22.81
## m0[16] 11.18
## m0[17] 22.49
## m0[18] 18.66
## m0[19] 17.31
## m0[20] 15.28
## y0[1] 6.71
## y0[2] 16.07
## y0[3] 17.33
## y0[4] 27.00
## y0[5] 15.21
## y0[6] 19.38
## y0[7] 12.65
## y0[8] 20.67
## y0[9] 7.39
## y0[10] 8.21
## y0[11] 8.83
## y0[12] 16.14
## y0[13] 24.05
## y0[14] 8.64
## y0[15] 22.00
## y0[16] 12.15
## y0[17] 23.97
## y0[18] 18.96
## y0[19] 16.16
## y0[20] 16.52
## m1[1] 8.65
## m1[2] 10.25
## m1[3] 11.86
## m1[4] 13.46
## m1[5] 15.06
## m1[6] 16.67
## m1[7] 18.27
## m1[8] 19.87
## m1[9] 21.48
## m1[10] 23.08
## y1[1] 7.33
## y1[2] 11.79
## y1[3] 11.44
## y1[4] 17.09
## y1[5] 15.45
## y1[6] 18.33
## y1[7] 17.94
## y1[8] 16.61
## y1[9] 17.38
## y1[10] 18.59
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -34.04 -33.73 1.22 1.14 -36.46 -32.63 1.00 614 1031
## b0 6.71 6.64 2.01 1.92 3.27 10.14 1.02 303 337
## b1 1.81 1.81 0.36 0.35 1.21 2.41 1.01 358 422
## s 3.70 3.61 0.66 0.60 2.77 4.95 1.00 1238 1161
## m0[1] 12.14 12.14 1.12 1.08 10.26 14.02 1.01 344 596
## m0[2] 14.99 14.98 0.85 0.84 13.64 16.45 1.00 735 1353
## m0[3] 15.31 15.29 0.84 0.82 13.98 16.74 1.00 848 1353
## m0[4] 21.09 21.10 1.32 1.30 18.94 23.22 1.00 981 1270
## m0[5] 16.84 16.83 0.85 0.84 15.49 18.29 1.00 2109 1587
## m0[6] 18.89 18.89 1.01 1.01 17.23 20.58 1.00 2030 1492
## m0[7] 14.15 14.13 0.90 0.85 12.73 15.68 1.00 512 994
## m0[8] 18.96 18.96 1.02 1.02 17.28 20.67 1.00 1962 1492
## m0[9] 10.51 10.49 1.36 1.29 8.23 12.83 1.01 313 424
## m0[10] 13.71 13.69 0.94 0.89 12.19 15.30 1.01 447 783
## m0[11] 11.34 11.31 1.23 1.18 9.30 13.43 1.01 325 473
## m0[12] 13.42 13.41 0.97 0.91 11.84 15.06 1.01 416 730
## m0[13] 22.97 22.98 1.62 1.60 20.29 25.54 1.00 736 1035
## m0[14] 8.85 8.82 1.63 1.55 6.10 11.66 1.01 305 339
## m0[15] 22.71 22.72 1.58 1.55 20.10 25.21 1.00 757 1071
## m0[16] 11.32 11.29 1.24 1.18 9.27 13.42 1.01 324 473
## m0[17] 22.39 22.40 1.53 1.50 19.87 24.82 1.00 789 1095
## m0[18] 18.64 18.64 0.99 0.98 17.03 20.29 1.00 2163 1537
## m0[19] 17.33 17.31 0.87 0.87 15.92 18.80 1.00 2254 1569
## m0[20] 15.34 15.32 0.84 0.83 14.01 16.77 1.00 859 1352
## y0[1] 12.19 12.23 4.01 3.89 5.57 18.86 1.00 1559 1789
## y0[2] 15.03 15.11 3.85 3.66 8.82 21.25 1.00 1882 2015
## y0[3] 15.30 15.25 3.79 3.51 9.06 21.68 1.00 2017 2014
## y0[4] 21.26 21.34 3.88 3.83 14.86 27.34 1.00 1984 1933
## y0[5] 16.63 16.67 3.84 3.61 10.31 22.73 1.00 2036 1786
## y0[6] 18.92 18.96 3.97 3.81 12.38 25.22 1.00 1923 1862
## y0[7] 14.04 13.99 3.94 3.72 7.58 20.55 1.00 1773 1845
## y0[8] 18.92 19.00 3.99 3.78 12.08 25.46 1.00 1890 1911
## y0[9] 10.47 10.50 3.98 3.81 3.94 17.22 1.01 1168 1856
## y0[10] 13.70 13.57 3.89 3.86 7.40 20.02 1.00 1986 2014
## y0[11] 11.40 11.29 3.90 3.88 4.85 17.77 1.00 1223 1452
## y0[12] 13.38 13.29 3.88 3.83 7.34 19.65 1.00 1499 1812
## y0[13] 22.92 22.82 3.99 3.96 16.53 29.80 1.00 1832 1835
## y0[14] 8.84 8.80 4.07 3.88 2.19 15.52 1.00 1023 1971
## y0[15] 22.78 22.86 4.11 3.93 16.09 29.33 1.00 1621 1881
## y0[16] 11.27 11.29 4.00 3.86 4.80 18.09 1.00 1320 1693
## y0[17] 22.34 22.36 4.12 4.14 15.70 28.96 1.00 1556 1830
## y0[18] 18.72 18.72 3.86 3.67 12.51 25.17 1.00 2232 1933
## y0[19] 17.43 17.42 3.81 3.72 11.19 23.56 1.00 1900 1802
## y0[20] 15.36 15.33 3.88 3.71 9.37 21.64 1.00 1930 1808
## m1[1] 8.85 8.82 1.63 1.55 6.10 11.66 1.01 305 339
## m1[2] 10.42 10.40 1.37 1.30 8.11 12.77 1.01 312 420
## m1[3] 11.98 11.98 1.14 1.09 10.06 13.91 1.01 340 544
## m1[4] 13.55 13.54 0.95 0.89 12.01 15.17 1.01 430 744
## m1[5] 15.12 15.11 0.84 0.84 13.79 16.56 1.00 779 1316
## m1[6] 16.69 16.68 0.84 0.82 15.36 18.14 1.00 1926 1594
## m1[7] 18.26 18.26 0.95 0.93 16.73 19.85 1.00 2250 1629
## m1[8] 19.83 19.83 1.13 1.12 17.95 21.70 1.00 1377 1369
## m1[9] 21.40 21.40 1.37 1.35 19.14 23.59 1.00 921 1190
## m1[10] 22.97 22.98 1.62 1.60 20.29 25.54 1.00 736 1035
## y1[1] 8.84 8.82 4.04 3.85 2.26 15.57 1.00 1165 1562
## y1[2] 10.43 10.51 4.05 3.89 3.87 17.15 1.00 1254 1649
## y1[3] 12.02 12.02 3.89 3.80 5.54 18.29 1.00 1800 1835
## y1[4] 13.59 13.66 3.95 3.73 7.10 20.20 1.00 1856 1858
## y1[5] 15.08 15.07 3.88 3.75 8.62 21.33 1.00 1899 1763
## y1[6] 16.60 16.60 3.80 3.62 10.50 22.83 1.00 1790 1750
## y1[7] 18.26 18.20 4.02 4.02 11.92 24.52 1.00 2078 1825
## y1[8] 19.72 19.76 3.98 3.87 13.06 25.98 1.00 1986 1657
## y1[9] 21.46 21.53 4.03 3.91 14.70 28.01 1.00 1893 1856
## y1[10] 22.88 22.94 4.12 3.89 15.99 29.64 1.00 1827 1651
mdl=cmdstan_model('./ex10.stan')
fn(mdl,data)## Initial log joint probability = -154.178
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 17 -18.0071 0.00122729 0.000938315 0.9814 0.9814 27
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -18.01
## b0 4.28
## b1 2.22
## s 0.77
## m0[1] 10.95
## m0[2] 14.46
## m0[3] 14.85
## m0[4] 21.96
## m0[5] 16.73
## m0[6] 19.25
## m0[7] 13.43
## m0[8] 19.34
## m0[9] 8.95
## m0[10] 12.88
## m0[11] 9.98
## m0[12] 12.53
## m0[13] 24.27
## m0[14] 6.91
## m0[15] 23.94
## m0[16] 9.95
## m0[17] 23.56
## m0[18] 18.95
## m0[19] 17.33
## m0[20] 14.88
## y0[1] 10.62
## y0[2] 15.25
## y0[3] 15.28
## y0[4] 21.34
## y0[5] 16.40
## y0[6] 19.13
## y0[7] 13.34
## y0[8] 19.76
## y0[9] 6.52
## y0[10] 12.88
## y0[11] 9.29
## y0[12] 14.90
## y0[13] 24.43
## y0[14] 8.82
## y0[15] 24.94
## y0[16] 10.16
## y0[17] 23.90
## y0[18] 18.99
## y0[19] 18.60
## y0[20] 14.41
## m1[1] 6.91
## m1[2] 8.84
## m1[3] 10.77
## m1[4] 12.69
## m1[5] 14.62
## m1[6] 16.55
## m1[7] 18.48
## m1[8] 20.41
## m1[9] 22.34
## m1[10] 24.27
## y1[1] 6.00
## y1[2] 8.16
## y1[3] 11.67
## y1[4] 17.05
## y1[5] 14.56
## y1[6] 16.14
## y1[7] 19.30
## y1[8] 19.88
## y1[9] 23.76
## y1[10] 23.95
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -19.86 -19.45 1.40 1.04 -22.71 -18.41 1.00 641 821
## b0 4.37 4.37 0.81 0.74 3.11 5.71 1.01 364 392
## b1 2.16 2.18 0.16 0.15 1.89 2.41 1.01 422 629
## s 0.98 0.94 0.31 0.28 0.57 1.53 1.00 852 748
## m0[1] 10.86 10.87 0.44 0.39 10.13 11.54 1.01 496 625
## m0[2] 14.27 14.30 0.36 0.35 13.63 14.81 1.01 1103 1116
## m0[3] 14.65 14.68 0.36 0.36 14.01 15.19 1.01 1249 1234
## m0[4] 21.56 21.65 0.65 0.64 20.38 22.49 1.00 1064 1208
## m0[5] 16.48 16.52 0.40 0.40 15.75 17.06 1.01 1968 1422
## m0[6] 18.93 19.00 0.50 0.50 18.02 19.66 1.00 1496 1405
## m0[7] 13.27 13.29 0.37 0.35 12.63 13.84 1.01 816 1030
## m0[8] 19.01 19.08 0.51 0.50 18.10 19.75 1.00 1472 1384
## m0[9] 8.91 8.92 0.53 0.47 8.06 9.75 1.01 426 520
## m0[10] 12.74 12.75 0.38 0.36 12.09 13.33 1.01 703 864
## m0[11] 9.91 9.92 0.48 0.43 9.12 10.66 1.01 453 547
## m0[12] 12.40 12.41 0.39 0.36 11.74 13.00 1.01 646 803
## m0[13] 23.80 23.90 0.80 0.78 22.37 24.92 1.00 871 1133
## m0[14] 6.93 6.93 0.65 0.58 5.89 7.96 1.01 390 415
## m0[15] 23.49 23.59 0.78 0.75 22.10 24.58 1.00 899 1162
## m0[16] 9.88 9.89 0.48 0.44 9.09 10.63 1.01 452 547
## m0[17] 23.11 23.21 0.75 0.73 21.76 24.17 1.00 936 1162
## m0[18] 18.63 18.70 0.49 0.49 17.76 19.34 1.00 1579 1354
## m0[19] 17.06 17.11 0.42 0.42 16.30 17.67 1.00 1974 1423
## m0[20] 14.68 14.71 0.36 0.35 14.04 15.22 1.01 1265 1234
## y0[1] 11.06 10.90 22.32 1.51 5.30 17.29 1.00 1777 1948
## y0[2] 12.12 14.27 59.72 1.50 8.01 19.74 1.00 1801 2059
## y0[3] 12.71 14.66 51.74 1.53 8.79 21.04 1.00 1995 1822
## y0[4] 22.08 21.59 31.86 1.63 15.55 28.10 1.00 1785 1957
## y0[5] 19.41 16.59 87.47 1.49 11.40 23.12 1.00 1872 1785
## y0[6] 19.16 18.98 31.37 1.60 12.45 25.25 1.00 1956 1851
## y0[7] 12.99 13.30 13.93 1.52 7.26 19.50 1.00 1845 1887
## y0[8] 12.29 19.04 291.05 1.60 12.60 25.92 1.00 2082 1664
## y0[9] 8.69 8.87 59.04 1.59 2.86 14.70 1.00 1754 1839
## y0[10] 12.12 12.73 25.18 1.51 5.66 18.23 1.00 1905 1891
## y0[11] 12.87 9.94 197.15 1.58 4.12 16.36 1.00 1891 1678
## y0[12] 13.57 12.39 76.38 1.46 4.81 18.67 1.00 2024 1794
## y0[13] 17.83 23.88 256.85 1.77 17.83 30.27 1.00 1738 1676
## y0[14] 6.22 6.97 52.92 1.69 0.51 14.66 1.00 1593 1892
## y0[15] 23.11 23.57 22.65 1.70 18.29 29.30 1.00 1832 1880
## y0[16] 11.25 9.92 93.64 1.61 2.12 16.64 1.00 1757 2015
## y0[17] 19.23 23.17 166.73 1.77 17.31 29.32 1.00 1707 1885
## y0[18] 18.78 18.67 14.50 1.47 12.54 24.78 1.00 2094 1959
## y0[19] 15.70 17.11 67.69 1.50 10.71 23.64 1.00 2102 1982
## y0[20] 17.56 14.76 47.88 1.49 9.04 21.81 1.00 2029 2041
## m1[1] 6.93 6.93 0.65 0.58 5.89 7.96 1.01 390 415
## m1[2] 8.80 8.81 0.54 0.48 7.94 9.65 1.01 424 520
## m1[3] 10.68 10.69 0.45 0.40 9.93 11.38 1.01 486 614
## m1[4] 12.55 12.57 0.38 0.36 11.90 13.15 1.01 671 825
## m1[5] 14.43 14.46 0.36 0.35 13.79 14.97 1.01 1158 1144
## m1[6] 16.30 16.34 0.39 0.40 15.59 16.88 1.01 1937 1380
## m1[7] 18.18 18.24 0.47 0.47 17.33 18.83 1.00 1715 1340
## m1[8] 20.05 20.13 0.56 0.55 19.04 20.85 1.00 1254 1343
## m1[9] 21.93 22.02 0.68 0.66 20.70 22.87 1.00 1030 1163
## m1[10] 23.80 23.90 0.80 0.78 22.37 24.92 1.00 871 1133
## y1[1] 8.44 6.91 50.65 1.58 1.01 12.73 1.00 1778 1858
## y1[2] 8.41 8.77 27.55 1.54 2.75 14.16 1.00 1672 1894
## y1[3] 16.95 10.69 274.03 1.47 4.35 16.88 1.00 1935 1750
## y1[4] -19.11 12.54 1210.82 1.44 4.63 17.09 1.00 1838 1440
## y1[5] 14.83 14.46 24.95 1.49 7.85 21.66 1.00 1940 1973
## y1[6] 13.39 16.31 74.00 1.57 9.63 22.88 1.00 1748 1994
## y1[7] 18.47 18.25 31.24 1.61 11.43 24.74 1.00 1857 2033
## y1[8] 22.11 20.10 60.92 1.67 13.71 26.85 1.00 1985 2017
## y1[9] 21.26 21.96 33.49 1.61 16.26 28.30 1.00 1852 1427
## y1[10] 24.19 23.84 33.26 1.78 17.41 30.22 1.00 1892 2014
censored data
objective variable has NA
(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)ex11-0.stan
data{
int N0;
array[N0] vector[2] xy;
int N1;
vector[N1] x1;
}
parameters{
vector[2] m;
cov_matrix[2] s;
}
model{
target+=multi_normal_lpdf(xy | m, s);
x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
vector[2] xy1;
xy1=multi_normal_rng(m,s);
real cr;
cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)## [1] 0.899qplot(x,y)
L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)## [1] 0.611qplot(x0,y0)
mdl=cmdstan_model('./ex11-0.stan')
data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -11806.3
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 70 -102.836 0.00249015 0.00278684 1 1 84
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -102.84
## m[1] 4.48
## m[2] 19.81
## s[1,1] 9.14
## s[2,1] 36.64
## s[1,2] 36.64
## s[2,2] 169.72
## xy1[1] 10.96
## xy1[2] 42.13
## cr 0.93mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.3 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.5 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -97.54 -97.20 1.70 1.50 -100.78 -95.51 1.00 430 727
## m[1] 4.53 4.52 0.61 0.62 3.57 5.55 1.01 590 777
## m[2] 21.26 21.17 5.92 5.75 11.45 31.29 1.01 136 103
## s[1,1] 11.75 11.24 3.52 3.20 7.22 18.09 1.00 795 1027
## s[2,1] 42.02 40.74 24.01 21.39 3.44 80.93 1.02 159 230
## s[1,2] 42.02 40.74 24.01 21.39 3.44 80.93 1.02 159 230
## s[2,2] 226.14 190.22 151.62 128.70 51.37 509.45 1.02 144 195
## xy1[1] 4.58 4.64 3.55 3.49 -1.46 10.27 1.00 1790 1846
## xy1[2] 21.26 22.47 15.84 13.87 -6.84 44.80 1.00 923 578
## cr 0.81 0.91 0.27 0.07 0.15 0.97 1.01 199 251
xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])## [1] 0.764qplot(xy[,,1],xy[,,2])
objective variable is censored
y i=1-N, y~N(m,s)
actual ya i=1-Na
lower censored yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
upper censored yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)
cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))
P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)ex11-1.stan
data{
int N;
vector[N] x;
vector[N] y;
real L;
int Nl;
vector[Nl] xl;
real U;
int Nu;
vector[Nu] xu;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
for(i in 1:Nl)
target+=normal_lcdf(L | b0+b1*xl[i],s);
for(i in 1:Nu)
target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)
xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitionmdl=cmdstan_model('./ex8-0.stan')
data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -337.674
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 32 -13.2188 0.000355531 7.26668e-06 1 1 45
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.22
## b0 15.37
## b1 1.74
## s 2.63
## m0[1] 22.77
## m0[2] 20.67
## m0[3] 26.59
## m0[4] 23.85
## m0[5] 26.79
## m0[6] 17.83
## m0[7] 23.75
## m0[8] 27.31
## m0[9] 22.03
## y0[1] 23.68
## y0[2] 19.18
## y0[3] 24.92
## y0[4] 24.52
## y0[5] 28.71
## y0[6] 15.77
## y0[7] 25.93
## y0[8] 26.79
## y0[9] 19.70
## m1[1] 15.74
## m1[2] 17.40
## m1[3] 19.05
## m1[4] 20.70
## m1[5] 22.35
## m1[6] 24.01
## m1[7] 25.66
## m1[8] 27.31
## m1[9] 28.97
## m1[10] 30.62
## y1[1] 14.98
## y1[2] 17.36
## y1[3] 17.70
## y1[4] 21.34
## y1[5] 22.11
## y1[6] 27.50
## y1[7] 23.26
## y1[8] 26.99
## y1[9] 25.34
## y1[10] 33.60mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.23 -13.74 1.76 1.33 -17.83 -12.43 1.01 284 323
## b0 15.03 15.07 4.58 3.44 8.18 22.24 1.02 264 218
## b1 1.81 1.80 0.90 0.69 0.38 3.25 1.02 275 246
## s 3.79 3.40 1.55 1.11 2.22 6.65 1.00 584 676
## m0[1] 22.71 22.74 1.45 1.19 20.37 24.88 1.01 606 658
## m0[2] 20.54 20.59 2.11 1.65 17.23 23.90 1.01 319 287
## m0[3] 26.68 26.62 1.98 1.63 23.66 30.09 1.01 642 740
## m0[4] 23.84 23.87 1.35 1.15 21.72 25.97 1.01 1280 979
## m0[5] 26.89 26.81 2.05 1.70 23.76 30.41 1.01 604 654
## m0[6] 17.58 17.66 3.38 2.60 12.45 22.92 1.02 272 220
## m0[7] 23.73 23.76 1.35 1.15 21.61 25.83 1.01 1205 950
## m0[8] 27.43 27.33 2.26 1.88 24.04 31.28 1.01 529 575
## m0[9] 21.95 21.99 1.63 1.32 19.33 24.47 1.01 430 389
## y0[1] 22.73 22.69 4.26 3.71 15.96 29.39 1.00 1839 1517
## y0[2] 20.59 20.79 4.62 3.97 13.23 27.60 1.00 925 838
## y0[3] 26.75 26.82 4.39 3.75 19.65 33.76 1.00 1369 1430
## y0[4] 23.75 23.81 4.37 3.63 16.76 30.47 1.00 1828 1773
## y0[5] 26.81 26.82 4.58 3.84 19.41 34.00 1.00 1241 1154
## y0[6] 17.56 17.49 5.05 4.23 9.44 25.78 1.00 562 756
## y0[7] 23.61 23.58 4.34 3.58 16.75 30.66 1.00 1633 1394
## y0[8] 27.38 27.45 4.59 4.11 19.93 34.44 1.00 1418 1393
## y0[9] 21.94 21.87 4.36 3.70 15.00 28.81 1.00 1475 1328
## m1[1] 15.42 15.46 4.39 3.34 8.89 22.27 1.02 265 223
## m1[2] 17.14 17.18 3.59 2.74 11.75 22.77 1.02 270 219
## m1[3] 18.85 18.91 2.81 2.15 14.54 23.34 1.02 282 235
## m1[4] 20.57 20.62 2.10 1.64 17.28 23.91 1.01 320 294
## m1[5] 22.28 22.32 1.54 1.26 19.78 24.67 1.01 489 491
## m1[6] 24.00 24.02 1.35 1.16 21.86 26.18 1.00 1393 1055
## m1[7] 25.72 25.70 1.66 1.38 23.08 28.53 1.01 966 809
## m1[8] 27.43 27.33 2.27 1.88 24.04 31.28 1.01 528 575
## m1[9] 29.15 29.02 3.00 2.49 24.60 34.16 1.01 409 451
## m1[10] 30.86 30.69 3.78 3.06 25.09 37.15 1.01 359 380
## y1[1] 15.39 15.37 5.94 4.74 6.43 24.41 1.01 431 491
## y1[2] 17.20 17.32 5.47 4.48 8.37 26.16 1.01 490 564
## y1[3] 18.90 18.93 4.88 4.10 11.40 26.59 1.00 695 758
## y1[4] 20.54 20.71 4.62 3.79 13.15 27.37 1.00 1068 1004
## y1[5] 22.37 22.48 4.32 3.63 15.38 29.01 1.00 1586 1601
## y1[6] 23.88 23.88 4.38 3.64 17.23 30.76 1.00 2033 1817
## y1[7] 25.59 25.67 4.36 3.50 18.77 32.36 1.00 1969 1722
## y1[8] 27.60 27.53 4.83 4.04 20.09 35.41 1.00 1270 1108
## y1[9] 29.28 29.13 4.89 4.13 21.78 37.01 1.00 984 1277
## y1[10] 30.78 30.82 5.63 4.62 22.27 39.42 1.01 641 798
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
data=list(N=n,x=xya$x,y=xya$y,
L=L,Nl=nl,xl=xyl$x,
U=U,Nu=nu,xu=xyu$x,
N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan')
mle=mdl$optimize(data=data)## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Initial log joint probability = -166.149
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## 24 -22.7103 0.00150634 0.000243648 0.9838 0.9838 33
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -22.71
## b0 8.96
## b1 3.28
## s 4.43
## m0[1] 22.89
## m0[2] 18.95
## m0[3] 30.08
## m0[4] 24.93
## m0[5] 30.46
## m0[6] 13.59
## m0[7] 24.74
## m0[8] 31.44
## m0[9] 21.50
## y0[1] 25.16
## y0[2] 17.33
## y0[3] 36.80
## y0[4] 19.44
## y0[5] 30.77
## y0[6] 15.86
## y0[7] 30.83
## y0[8] 32.83
## y0[9] 21.00
## m1[1] 9.67
## m1[2] 12.78
## m1[3] 15.89
## m1[4] 19.00
## m1[5] 22.11
## m1[6] 25.22
## m1[7] 28.33
## m1[8] 31.45
## m1[9] 34.56
## m1[10] 37.67
## y1[1] 15.56
## y1[2] 21.52
## y1[3] 16.84
## y1[4] 18.80
## y1[5] 21.75
## y1[6] 22.76
## y1[7] 29.49
## y1[8] 30.55
## y1[9] 29.59
## y1[10] 28.64mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,'m',ch=T)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -23.11 -22.71 1.55 1.23 -26.02 -21.44 1.02 287 481
## b0 5.84 6.63 5.40 4.43 -3.62 12.30 1.02 220 281
## b1 3.99 3.80 1.14 0.93 2.63 5.92 1.01 240 255
## s 6.51 5.94 2.66 1.95 3.71 11.17 1.01 405 602
## m0[1] 22.81 22.91 1.91 1.70 19.64 25.65 1.00 751 670
## m0[2] 18.01 18.36 2.48 2.11 13.64 21.21 1.01 314 413
## m0[3] 31.57 31.15 2.97 2.50 27.72 36.81 1.01 641 526
## m0[4] 25.29 25.24 1.96 1.76 22.37 28.45 1.00 1344 882
## m0[5] 32.03 31.61 3.07 2.58 28.06 37.49 1.01 604 504
## m0[6] 11.48 12.07 3.94 3.27 4.41 16.17 1.01 229 287
## m0[7] 25.06 25.02 1.95 1.76 22.14 28.15 1.00 1290 938
## m0[8] 33.22 32.72 3.34 2.79 28.97 39.27 1.01 529 490
## m0[9] 21.11 21.30 2.03 1.77 17.67 23.96 1.01 508 523
## y0[1] 23.03 22.93 7.07 5.98 12.04 34.19 1.00 1770 1714
## y0[2] 17.97 18.27 7.47 6.37 5.61 28.79 1.00 1334 973
## y0[3] 31.35 30.98 7.66 6.55 19.90 43.92 1.01 1782 1284
## y0[4] 24.99 24.94 7.09 6.20 14.32 36.51 1.00 1983 1576
## y0[5] 32.06 31.63 7.45 6.10 20.99 44.27 1.00 1701 1715
## y0[6] 11.19 12.10 8.06 6.68 -2.74 21.93 1.00 658 965
## y0[7] 25.03 25.01 7.28 6.18 13.60 36.42 1.00 2005 1688
## y0[8] 33.22 32.78 7.88 6.66 21.63 46.19 1.00 1605 968
## y0[9] 21.22 21.41 7.37 6.01 9.65 32.51 1.00 1833 1680
## m1[1] 6.71 7.48 5.17 4.25 -2.33 12.86 1.02 221 276
## m1[2] 10.50 11.12 4.18 3.49 3.12 15.44 1.02 226 294
## m1[3] 14.29 14.79 3.26 2.72 8.32 18.29 1.01 244 324
## m1[4] 18.07 18.42 2.46 2.11 13.73 21.27 1.01 316 418
## m1[5] 21.86 22.02 1.96 1.71 18.57 24.69 1.01 598 564
## m1[6] 25.65 25.59 1.99 1.79 22.70 28.88 1.00 1413 857
## m1[7] 29.44 29.12 2.53 2.23 26.05 33.79 1.01 940 733
## m1[8] 33.23 32.72 3.34 2.79 28.98 39.28 1.01 528 490
## m1[9] 37.02 36.29 4.28 3.55 31.68 44.70 1.01 404 450
## m1[10] 40.81 39.91 5.26 4.25 34.31 50.17 1.01 351 397
## y1[1] 6.82 7.81 8.84 7.44 -8.51 18.62 1.01 596 590
## y1[2] 10.33 10.85 8.36 6.66 -2.84 21.60 1.01 692 685
## y1[3] 14.57 14.98 7.78 6.79 1.78 26.13 1.00 1172 1116
## y1[4] 18.06 18.30 7.67 6.07 6.18 29.56 1.00 1775 1200
## y1[5] 22.03 22.35 7.32 6.40 10.10 33.31 1.00 1714 1227
## y1[6] 25.58 25.49 7.19 6.04 14.84 37.09 1.00 1947 1693
## y1[7] 29.23 29.03 7.32 6.27 18.33 41.06 1.00 2056 1366
## y1[8] 33.18 32.64 7.47 6.40 21.99 46.14 1.00 1692 1351
## y1[9] 37.24 36.42 8.56 6.99 24.87 51.54 1.01 1213 1103
## y1[10] 40.54 39.77 8.42 7.04 28.84 54.86 1.00 785 988
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
sensitivity/specificity analysis
sensitivity: true positive rate TPR = TP/(TP+FN)
specificity: true negative rate TFR = TN/(FP+TN)
ROC curve: se vs 1-sp
positive predictive value ppv = TP/(TP+FP)
negative predictive value npv = TN/(TN+FN)ex14.stan
estimating sens and spec
data {
int N;
array[N] int x;
array[N] int y;
}
parameters {
real<lower=0,upper=1> p;
real<lower=0,upper=1> se;
real<lower=0,upper=1> sp;
}
model {
p~uniform(0,1);
se~uniform(0,1);
sp~uniform(0,1);
for (i in 1:N) {
y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
}
}
generated quantities {
real ppv;
real npv;
ppv=se*p/((se*p)+((1-p)*(1-sp)));
npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}n=20
x=sample(0:1,n,replace=T)
p=(x+0.5)*0.5
y=rbinom(n,1,p)
data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex14.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -14.7148
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -13.4428 0.000615617 0.000249883 1 1 9
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.44
## p 0.88
## se 0.58
## sp 0.63
## ppv 0.92
## npv 0.17mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -19.30 -18.95 1.32 1.03 -21.77 -17.89 1.00 755 840
## p 0.49 0.48 0.28 0.36 0.06 0.94 1.01 710 423
## se 0.57 0.58 0.13 0.14 0.35 0.77 1.00 3177 1523
## sp 0.59 0.61 0.15 0.16 0.34 0.82 1.00 1820 1245
## ppv 0.55 0.58 0.29 0.37 0.08 0.96 1.01 742 556
## npv 0.57 0.60 0.28 0.36 0.07 0.96 1.01 749 594
ppv=mcmc$draws('ppv')
npv=mcmc$draws('npv')
qplot(ppv,npv)
2x2 cross table
Effect occur y=1 with probabilty p01, p11 from each cause x{0,1}
event frequncy nxy of effect y{0,1} by cause x{0,1}
n01~B(n0.,p0)
n11~B(n1.,p1)
n01=n0p0, n00=n0(1-p0)
n11=n1p1, n10=n1(1-p1)
p00=n00/n=n0(1-p0)/n, p01=n01/n=n0p0/n
p10=n10/n=n1(1-p1)/n, p11=n11/n=n1p1/n
Cramer'V (chi2/n/(min(row,column)-1))^.5
in 2x2
crv =(n11*n00-n10*n01)/(n0.*n1.*n.0*n.1)^.5
=(n0(1-p0)n1p1-n0p0n1(1-p1))/(n0n1(n0(1-p0)+n1(1-p1))(n0p0n1p1))^.5
kappa coefficient k=(po-pe)/(1-pe)
po: Observed agreement (proportion of times both raters agreed)
pe: Expected agreement under independence
po=p00+p11
=(n0(1-p0)+n1p1)/n
pe=(p00+p01)(p00+p10)(p10+p11)(p01+p11)
=n0/n*(n0(1-p0)+n1(1-p1))/n*(n0p0+n1p1)/n*n1/nex16-1.stan
data {
int n;
int n0;
int n01;
int n1;
int n11;
}
parameters {
real<lower=0,upper=1> p0;
real<lower=0,upper=1> p1;
}
model {
n01~binomial(n0,p0);
n11~binomial(n1,p1);
}
generated quantities {
real RR;
RR=p1/p0;
real OR;
OR=(p1/(1-p1))/(p0/(1-p0));
}ex16-2.stan
data {
int n;
int n0;
int n01;
int n1;
int n11;
}
parameters {
real<lower=0,upper=1> p0;
real<lower=0,upper=1> p1;
}
model {
n01~binomial(n0,p0);
n11~binomial(n1,p1);
}
generated quantities {
real RR;
RR=p1/p0;
real OR;
OR=(p1/(1-p1))/(p0/(1-p0));
real crv;
crv=(n0*(1-p0)*n1*p1-n0*p0*n1*(1-p1))/(n0*n1*(n0*(1-p0)+n1*(1-p1))*(n0*p0+n1*p1))^.5;
real k;
real po;
po=(n0*(1-p0)+n1*p1)/n;
real pe;
pe=n0/n*(n0*(1-p0)+n1*(1-p1))/n*(n0*p0+n1*p1)/n*n1/n;
k=(po-pe)/(1-pe);
}n0=30
n01=rbinom(1,n0,0.3)
n1=30
n11=rbinom(1,n1,0.6)
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)
mdl=cmdstan_model('./ex16-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -55.3333
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -40.7173 0.00355404 6.49073e-05 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -40.72
## p0 0.40
## p1 0.43
## RR 1.08
## OR 1.15mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -44.48 -44.20 0.89 0.66 -46.25 -43.58 1.00 939 972
## p0 0.40 0.40 0.08 0.08 0.27 0.54 1.00 2047 1583
## p1 0.44 0.44 0.08 0.09 0.30 0.58 1.00 2165 1521
## RR 1.13 1.08 0.33 0.31 0.67 1.73 1.00 2068 1567
## OR 1.30 1.14 0.67 0.55 0.51 2.52 1.00 2109 1556
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)
mdl=cmdstan_model('./ex16-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -58.2846
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 4 -40.7173 0.0154814 0.000587882 0.9942 0.9942 7
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -40.72
## p0 0.40
## p1 0.43
## RR 1.08
## OR 1.15
## crv 0.03
## k 0.49
## po 0.52
## pe 0.06mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -44.48 -44.20 0.89 0.66 -46.25 -43.58 1.00 939 972
## p0 0.40 0.40 0.08 0.08 0.27 0.54 1.00 2047 1583
## p1 0.44 0.44 0.08 0.09 0.30 0.58 1.00 2165 1521
## RR 1.13 1.08 0.33 0.31 0.67 1.73 1.00 2068 1567
## OR 1.30 1.14 0.67 0.55 0.51 2.52 1.00 2109 1556
## crv 0.03 0.03 0.12 0.12 -0.17 0.22 1.00 2111 1504
## k 0.49 0.48 0.06 0.06 0.38 0.58 1.00 2111 1553
## po 0.52 0.52 0.06 0.06 0.42 0.61 1.00 2115 1628
## pe 0.06 0.06 0.00 0.00 0.06 0.06 1.00 2065 1609
point of subjective equality PSE
PSE: 50% threshold for sensing the difference between two stimuli is equal
JND: Just noticeable difference, difference between 50% threshold and 75%
r~B(n,p) #reaction for stimuli
p=1/(1+exp(-(a+b*x)))
x=x1-x0, x0,x1 is strength of stimuli
PSE=-a/b
JND=(log(0.75/0.25)-a)/b-PSEex6-3-0.stan
mulit logistic regression
data{
int N;
int m;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
}
model{
y~binomial_logit(m,b0+b1*x);
}n=20
m=10
x=runif(n,-2,2)
y=rbinom(n,m,1/(1+exp(-(-2+3*x))))
glm(y/m~x,family=binomial('logit'))##
## Call: glm(formula = y/m ~ x, family = binomial("logit"))
##
## Coefficients:
## (Intercept) x
## -2.32 2.78
##
## Degrees of Freedom: 19 Total (i.e. Null); 18 Residual
## Null Deviance: 16.6
## Residual Deviance: 1.9 AIC: 10data=list(N=n,m=m,x=x,y=y)
mdl=cmdstan_model('./ex6-3-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -286.995
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 10 -52.0684 0.000462636 0.00031605 1 1 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -52.07
## b0 -2.32
## b1 2.78mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -53.07 -52.77 0.98 0.71 -55.09 -52.12 1.00 751 791
## b0 -2.44 -2.42 0.44 0.46 -3.22 -1.76 1.00 670 688
## b1 2.92 2.88 0.43 0.43 2.29 3.71 1.00 660 722
b0=mcmc$draws('b0') |> as.vector()
b1=mcmc$draws('b1') |> as.vector()
pse=-b0/b1
quantile(pse,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)## 0% 2.5% 5% 25% 50% 75% 95% 97.5% 100%
## 0.531 0.665 0.690 0.778 0.837 0.896 0.981 1.009 1.196jnd=(log(0.75/0.25)-b0)/b1-pse
quantile(jnd,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)## 0% 2.5% 5% 25% 50% 75% 95% 97.5% 100%
## 0.247 0.287 0.296 0.344 0.382 0.420 0.480 0.508 0.694x1=runif(length(b0),-2,2)
p=1/(1+exp(-(b0+b1*x1)))
pm=1/(1+exp(-(median(b0)+median(b1)*x1)))
qplot(x1,pm,col=I('darkgray'),ylab='p')+
geom_line(aes(x=x1,p=p),col='red')